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INTRODUCTION 



I . 



A. BACKGROUND 

The Marcum- Swerling (M-S) models are the most commonly 
used target models in modern radar analysis. The other target 
models have also developed over period of time, such as chi- 
square (Weinstock) , log-normal. However in this thesis only M- 
S models will be discussed. 

Histrocially , the earliest descriptions of a target were 
in terms of a single cross-sectional area value. This quantity 
was usually some type of an average cross section over the 
aspect angles which the system designer considered most 
probable. The approach led the system designers to associate 
a single cross-sectional area of unvarying value, with a 
target and thus led to the model of a steady-state target. 
Following this concept, systems were designed to achieve a 
somewhat arbitrarily specified signal-to-noise ratio at some 
specified maximum operating range. The inherent basic 
assumption was that the required detection capability and 
parameter estimation accuracy could be achieved if this goal 
was met. The dominant rule of the human observer in early 
detection systems tended to make this approach acceptable. As 
the requirements on radar systems became more demanding, and 
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as automatic processing of radar returns was developed, the 
situation underwent a change. The desire to optimize radar 
designs established the need for more precise target models. 

The initial work of J. I. Marcum of the RAND Corporation in 
the late 1940s applied the work of Rice of the Bell Telephone 
Laboratories, relating to steady-state signals immersed in 
noise, to the problem of radar signal detection. Marcum 
effectively gave a complete treatment to the statistical 
problem of a group of constant amplitude signal pulses in the 
presence of noise. His work resulted in an evaluation of 
previously untabulated functions and a direct application of 
the results to a gamut of problems in automatic detection. 
The statistical approach used in Marcum's studies is the basis 
for much of the work that followed in the radar detection 
area . 

The detection of signals resulting from a fluctuating 
target is basically different than for signals resulting from 
a non-fluctuating target. A requirement remained, therefore, 
to produce the same type of analysis and supporting tabulation 
of basic functions for a fluctuating target as had been done 
by Marcum for the nonfluctuating target. This work was done by 
Swerling, in the early 1950s. He extended Marcum's approach by 
employing four target models and two different density 
functions in conjunction with two extremes of correlation. 

Swerling ' s case 1 and case 2 are the target models which 
describe large complex targets, such as aircraft, rain 
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clutter, and terrain clutter. They represent two extremes of 
correlation, and the statistical model is derivable from the 
physical scattering characteristics of such bodies. It can 
also be stated that empirically derived data on such targets 
are in very good agreement with those obtained from the 
mathematical model. 

At the time Swerling was doing his work, it was realized 
that the model suitable for large complex targets would not 
give an adequate description of large simple structures. The 
problem of selecting a model to describe this type or class of 
target is a difficult one to handle directly. There is no 
obvious parallel development for it from the scattering 
characteristics of a specific body type. The approach 
employed was to establish qualitatively the basic statistical 
behavior of the target cross sections of interest. Having done 
this, a convenient form of distribution was selected. 

Swerling picked a class of distributions know as the chi- 
square class, which has the exponential density function 
distribution as one member. By an appropriate selection of a 
class parameter, namely, the number of degrees of freedom, the 
qualitative properties desired for the distribution associated 
with large simple structures was obtained. Swerling ‘s cases 3 
and 4 represent targets that behave as if they have four 
degrees of freedom and are valid for targets such as rockets, 
missiles, and space-based satellites. 

Swerling' s two classes of density functions evaluated at 
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two extremes of correlation, together with Marcum's constant 
target model (case 5), have been the bases of virtually all 
radar detection analyses. 

Five target models according to Marcum-Swerling scheme are 
as follows: 

a) Swerling case 1 - The echo pulses received from a 
target on any one scan are of constant amplitude throughout 
the entire scan but are independent ( uncorrelated ) from scan 
to scan. This assumption ignores the effect of the antenna 
beam shape on the echo amplitude. The probability density 
function of the radar cross section a is given by the 
exponential density function. 

w(o,o)=ie , " < ' /5| ( o>0 (1) 
o 

where <? is the average radar cross section 

b) Swerling case 2 - With the same density function as 
case 1 but the f luctulations are more rapid than in case 1 and 
are taken to be independent from pulse to pulse rather than 
from scan to scan. 

c) Swerling case 3 - The fluctuation is assumed to be 
independent from scan to scan, as in case 1, but the 
probability density function is given by the chi-square 
distribution with four degrees of freedom. 

w(o,o) = , 4 (Ko/o)*- 1 e~ 2a/ °, (K= 2) (2) 
V-K“l; • a q 2 



4 



d) Swerling case 4 - With the same density function but 
the fluctuation is pulse- to-pulse according to Eq(2). 

e) Marcum's model (case 5) - With the constant density 

function which represents steady-state target. 

The probability density function of equation (1), applies 
to a complex target consisting of many independent scatterers 
of approximately equal echo areas. The probability density 
function assumed in case 3 and 4 is more indicative of targets 
that can be represented as one large reflector together with 
other small reflectors. 

There are curves available that can be used to calculate 
the probability of detection for each of Swerling cases, but 
for parameters between those charted, the designer has to 
interpolate. This can be inaccurate sometimes. Therefore it is 
convenient and useful to provide accurate programs to 
calculate for any parameters needed by the user. 

B. RADAR RECEIVER MODEL 

1 . DESCRIPTION 

The typical super heterodyne radar receiver and the 
mathematically receiver model is depicted in Fig 1. The 
difference between square-law and linear envelope detector is 
that a square-law envelope detector is used in the small 
signal optimum receiver and a linear envelope detector is used 
in the large-signal receiver. For mathematical convenience, 
the square-law detector is applied, but the performance of 
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both detectors is practically the same for all signal- to-noise 
ratios . 
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Figure 1: Radar Receiver Model 

The multiple-pulse detection is used to improve the 
detectability of the signal by achieving integration gain. The 
noncoherent integration provides an integration gain even when 
the signal has a random phase and is rapidly fluctuating, such 
as Swerling case 2 and 4. By contrast, coherent integration 
needs a nonfluctuating or slowly fluctuating signal with 
predictable phase characteristics. 
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2 . ASSUMPTIONS 



The assumptions embodied in the M-S detection problems are 
as follows : 

• A received pulse train consisting of N samples of noise 
or signal-plus-noise is available. 

• The signal is imbedded in white Gaussian noise of known 
spectral density of N 0 /2 . 

• The signal is of unknown phase type, where the RF phase 
between pulses in the train are randomly distributed. 

• The processing consists of a matched filter and square- 
law envelope detector which operates on each pulse of the 
train and a linear integration which combines the n square-law 
envelope detected pulse. 

• The directivity pattern of the antenna is rectangular so 
that the n return pulses resulting during the antenna's dwell 
time on the target are unaffected by the antenna's radiation 
pattern . 

• The target's cross section can be described by a chi- 
squared distribution with 2k degree of freedom. 



w(o, o) 



1 K <Ko, K -' ^-(Kcfo) 

TK^iTJ oo 



(3) 



C. RADAR DETECTION PHILOSOPHY 

Radar detection is complicated by the fact that the target 
cross section a is a random variable fluctuating with time and 
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both noise and clutter. Extraction, or parameter estimation, 
is likewise complicated by random fluctuations of the target 
echo. For the radar detection case Neyman- Pearson criterion is 
used which needs neither a priori probabilities nor cost 
estimates. In radar terminology whose objective is to maximize 
the probability of detection for a given probability of false 
alarm. This objective can be accomplished by using a 
likelihood ratio test. Specially, there exists some 
nonnegative number T| such that hypothesis H : (i.e., target is 
present) is chosen when 

A( y) = * n <*> 

and hypothesis H 0 (no target present) is chosen otherwise. 
There are two types of errors which can be made. If we say a 
target is present when in fact it is not, an error of the 
first kind is made (see Fig 2). That is, we choose H : given 
that H 0 is true which is denoted as the probability of false 
alarm. Similarly an error of the second kind is made when H 0 
is chosen and in fact H x is true which is probability of miss. 
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To optimally detect the signal with unknown phase, 
consider the narrowband signal of duration T given by 



s ( t) =Aa ( t) cos [w 0 1+0 ( t) +<J>] 
=s I costy -s^int}) 



(5) 



where S z =Aa ( t) cos [w o t+0 ( t) ] , S q =Aa ( t ) sin [ w o t+0 ( t ) ] , A is signal 
amplitude, a(t) is an envelope function of duration T, 0(t) is 
the signal phase modulation. Waveform r(t) is assumed to be 
prefiltered by an ideal low-pass filter that is distortionless 
within (~w c ,w c ) and zero outside the interval; therefore the 
filter will not affect the band-limited input signal s(t) but 
results in a band-limited statistical process. Under this 
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assumption the input waveform r(t) can be described by the 
sampling theorem as 

r i ( t) = s i (t) + J2 i (t) i = 1,2, 3,. ..,2 f c T (6) 

Where 2 f c T is the number of samples, and t=i/2/ c . 

The hypothesis testing can be described as: 

H, : r = 8 + n ( target present ) 

(7) 

H 0 : r = n ( target absent ) 

Then under each hypothesis the probability density functions 
can be written as: 



£ 



U rSi )‘ 



f„+n(r) =- 



2oZ 



(v/2Wo n ) 



f n (r)=- 



S' 



20 : 



(8) 



(v/2«o n ) 



The likelihood ratio A(r) is given by 



A(r) = 









lim 



2f c T r , 2 

/ l 

eX P { "L oV A r } 

ft 2f C N 0 



f n (r) 

exp- 



2f c T 

{_ — 2 \ 
l-l 



2^o 

exp{--i- f T [r( t) -s( t) ] 2 dt) 

Nr, JO 



exp f T r(t) 2 dt} 

N 0 Jo 



(9) 



= exp [~-^-[ T s 2 (t) dt + -£- f T r(t) s(t) dt] 

N 0 J o N 0 Jo 

where E = J 0 T s 2 (t)dt is the energy of the signal. Substituting 
Eq(5) into Eq(9) yield the following expression for the 



10 



likelihood ratio A(r | <J) ) 



A (r|4>) = exp{--^+-J- [ f T r ( t) s i ( t) dt] cos4> 



¥ 0 ^0 



( 10 ) 



- -J- [f T r(t) s (t) dt] sin(J> 
W 0 Jo g 



Let the signal- to-noise ratio <J) -E/N 0 and 



y z (T)= — — f T r(C) Sj(C) dt 
M 4 /ifr Jo 



V'P 0 



y 0 (T) = — — f r(t)s p (t) dt 
^oV^ Jo 

using Eq(ll), Eq(10) can be written as follows 



(ID 



A(r|<j>) = e * /2 exp {\f^[y I [T) cos<J> - y 0 (T) sin<t>] } 



= e + /2 exp {y/ijr i (T) cos (<J>+a) } 



( 12 ) 



where a=tan _1 (y Q /y t ) , r( T) = [y Q ( T) +y z ( T) ] 1/2 is the envelope of the 
radar signal out of a filter matched to the waveform (2/N 0 
\\/ 1/2 )s(t), sampled at time t. Assume <{> is a uniformly 
distributed density function U(0 2k) . Averaging with respect 
to <|> yields the averaging likelihood ratio as follows 



A(r) =e* /2 — f 2n e ^ r{T)coe{ ^* a) d<J) 
2k Jo 



(13) 



=e~* /2 I 0 [vAJrr ( T) ] 



where I 0 (\y 1/2 r (T) ) is the modified Bessel function of the first 



11 



kind and order zero. Therefore the threshold test becomes : 



"i 

J 0 W r(T))] < e»' 2 A(r) (14) 

"o 

where e^A (r) is the operating threshold. The corresponding 
receiver implementation is shown in Fig 3 




Signal Signal 

p r aa m abaant 



Figure 3 : Optimum detector for radar signal of unknown phase 
D. OBJECTIVE 

The Marcum- Swer ling target model is the most commonly 
used model to calculate radar range performance. The objective 
of this thesis is to review the underlying theory of radar 
detection for this model and then develop a MATLAB programs to 
compute probability of detection and maximum detection range. 
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II. MODELS AND MATHEMATICAL METHODS 



The procedure for the calculation of average probability 
of detection P D for Marcum and Swerling cases is as follows : 

1) Find the single pulse characteristic function, C^p), by 
transforming the single pulse probability density function 
f s+n ( y ) • 

2) Find the N pulse characteristic function, C N (p)=[C 1 (p) ] N . 

3) Average the N pulse characteristic function over the target 
distribution to find the average of C N (p) ( C N (p) ) . 

4) Transform C N (p) to find the average ensemble probability 
density function f s . n ( y) and f n (y) for the impulse return. 

5) Find probability of false alarm P fa by integrating f n ( y) 
from y b to infinity. 

6) If P fa is given find y b by means of the mathematical 
recursive method. 

7) Find the average probability of detection P d , by integrating 
the density function of signal plus noise £ s +n(y) from the 
threshold y b to « . 

A. MARCUM'S (NON- FLUCTUATING) STEADY-STATE TARGET MODEL 

The non-fluctuating target model is applied to spherical 
or nearly spherical objects, such as balloons, many 
wavelengths in diameter. The target model, therefore, can be 
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represented by a constant-valued radar cross section. 

The Rayleigh and Rician probability density function of 
the envelope detected output for N pulses are given by 



f n^ r iW = r i e ~ r> ! r i*° 

- r £ e 2 I 0 (r i} ffi) ; r^O 



( 15 ) 



where \\i=E/N 0 is the single pulse peak signal- to-noise ratio 
for a steady target. Since the noise received in the ith 
observation interval is assumed to be statistically 
independent of the noise in every other observation interval, 
and the signals are noncoherent between the different 
observation interval, the initial phase (jh in Eq(5) in the ith 
interval is also statistically independent of <() j for all i^j . 
The likelihood ratio test in Eq (9) can be written as: 



T r ~\ - f i( r i> r 2' • • - r J s i' s 2 - • - s n) 

foTTTjo) 



_ £ f i( r A s i) 
' i-i fJFJo) 



= II e~* /2 I 0 (r iV /iJ7) 

i- 1 

= e-^ni 0 (r^) 

i-l 

From Eq(16), the test statistics is 



( 16 ) 



N 



H, 



lnJ 0 (r iV /iJF) < e* /2 A(r) 



2-1 



H n 



( 17 ) 
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For small signal, the modified Bessel function I 0 (x) can be 
approximated for x<l as: 



x 2 , x 4 



I Q (x) = i + — +— + 



64 



x 2 . x 4 



In T 0 (x) = ln(l + ^-+^- + . . . ) = (x 4 ) 

4 64 4 



(18) 



Therefore Eq(17) can be rearranged and modified as 



N 



rhi Y 



I* 1 



(19) 



where T|' is the operating threshold which is determined by 
specifying the false alarm probability. 

Marcum utilizes a square law detector law which allows the 
signal-plus-noise probability density function to be expressed 
directly in the signal-to-noise ratio y. To simplify the 
calculation, let 



AT 2 N 

y - E T - ' £ 

1*1 z i^l 



(20) 



which Y is compared to a suitably modified threshold Y h . 

To find the probability density of Y, first the 
probability density of f^tqi) is found by Jacobain 
transformation f s+n (r i ) from Eq(15) into f s + n (Q[i) as follows: 



<Ji*1 



f /_ \ _ _ 2 



I 0 (yJTq^) 



q^o 



( 21 ) 
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Since random variable Y is the sum of N statistically 
independent random variables q it by applying the relationship 
between the sum random variable and the individual random 
variable of the characteristic function, and the Campbell and 
Foster tables of Fourier transforms ; the C Y (p) is : 



C Y {p ) = H C (p) = 0 f e Jpg ‘f s , N (q i ) dq i 

i-1 1 



N r c 

= n / 

1 = 0 J 0 



,-Wf/2 



JPQi + 

W* 






( 22 ) 



(1+p) w 

From Campbell and Foster Tables of Fourier transforms pair 
650.0, the probability density function of f s+n (Y) can be found 
by taking the inverse transform of C Y (p) . 



= (■%£) iN ' 1)/2 e'™* /2 I^iy/mfY) , y^O (23) 

For small Y, the modified bessel function I N . 1 (Y) can be 
approached as : 

J m( y)=^L [1 + Zl + Zl ♦ ...] (24) 

2 m w\ 2 2 (m+1) 2 5 (/n+1) (w+2) 



Therefore the resulting normalized square-law detected 
probability density function is 






YN-i e -r-Mf/2 

(N 1 1) ! ' 



y^o 



( 25 ) 
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for noise only the probability of density function f n ( y) is 

yN- \~-Y 

f ” m ’ ITFTTi' Yi0 (36) 

A computer simulation of f n (Y) at the square law detector 
output is given in Figure 4 : 




Figure 4 : Simulation of noise at suare law detector output 

After square-law detection, the normalized square-law 
detected variate (y) is summed over the N-pulses in the pulse 
train and then compared against a normalized threshold voltage 
(y b ) to determine the presence or absence of a target return. 

The probability of false alarm can be obtained by 
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integrating from the threshold (y b ) to infinity 



Pfa = 



■/. 



yN-l e -y 

yb (N~ 1 ) ! 



dY 



( 27 ) 



The incomplete Pearson gamma function which is very useful 
for describing M-S model is defined as: 

Ku.pt ( 28 ) 

Jo pi 

And in term of the Pearson's incomplete gamma function, Eq(27) 
is 



P fa (N,Y b ) = 1-IlYjN 1 ' 2 , (N-l) ] ( 29 ) 

where u= (y b /N 1/2 ) and P = N-l . 

Eq(29) can be successively integrated by parts and given 

as : 



Pfa y b ) 



y*2£!l \ 1+ (N-l) + (N-l) (N- 2) + , 

(N-l) ! Y b Y b 2 



N-l 

-E 

K= 0 



Yb 



K e' Yb 



K\ 



( 30 ) 



The above equation also can be approximated for N>>1 by 
letting 



r x j. N - 1 x (N-l) (N- 2) + i _ 1 

Y t Y b 2 " 

and by applying Stirling's approximation 



N >1 



( 31 ) 



AT N 

Nl =SZnN(~) 
e 

The probability of false alarm in Eq(30) 



( 32 ) 

can be approximated 



18 



as 



E> <■» 

*fa 



e‘ Kfc+w 



N> 1 



(33) 



2n N Y-b-N+1 
The probability of detection for steady- target (case 5) 
can be obtain by integrating f s . n (Y) from y b to •*» and is given 
by 



P d (Y) = f“ f stn (Y) dY 

J Y b 

= i-J'J'ijl) w-D/2 e -H*»)/2j^ (v /5mJTP ) r;>o 



(34) 



The Q function is defined as 



0 N (cc,P)=fp W-^) exp [ - ( *1^1 ) ] i N _ 1 (av)dv (35) 

For N pulses the probability of detection P D can be 
written in terms of Q function as follows: 

P D (y) =Q{SZm >s [Tr b ) (36) 

Another approximation can be made by using Gram-Charlier 
series expansion ( Appendix A ) and noting that Gaussian 
family is closed under linear operation. 

m i = (-J) dC ^p P> lp-o * W(l++/2) • 

m 2 - (-j) ; d ^ r( „ P) | p .„ - W 2 (l^/2) 2 *W(l*t) <3?) 

dp N 

a 2 =m 2 -ml = N{l+\\i) 



where m i is the ith moment of the distribution of Y and c is 
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the variance of Y. 



f ( Y )s 

J -8+n v ' 



[r-W(l+t/ 2 ) ] : 



/SnN ( l _ +i[r )" 2iV( 1 +ijf ) 



f(Y) a- 1 e ( ~ r - M2/2W 



(38) 



Therefore the probability of false alarm can be represented 
as : 



fa 



* (— 

J * /O tt 



_JY-N If 






y b sfTizN 



* *[ Yb ~ N ' 



2W dy =< f — — e 2 dZ where Z= — — 
•SR 



(39) 



s[N 



Y b can be solved from the above equation 



Y b * y[R$-HP fa ) +N 

The probability of detection can be obtained as 



(40) 



. , (Y b -N (l+T/2)) 

Pd a 4>[ ; ■ 

v/Wd + T) 



] 



(41) 



B. SWERLING'S (FLUCTUATING) TARGET MODELS 
1. SWERLING CASES 1 AND 2 

The radar cross section of a target (a) is the area of a 
hypothetical reflector that scatters all radar beam energy it 
intercepts omnidirectionally, such that it produces an echo at 
the radar antenna equal to that energy actually received from 
the target; that is 
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a= (power received by antenna) 

( incident power density per 4 71 rad) 

( 42 ) 

=lim 4jii? 2 -^L 

m 2 

Where is the electric field incident on the target, E R is 
the reflected field as measured at the radar receiver antenna, 
and R is the distance from the target to the radar receiver 
antenna . 

The total electric field received from a complex target 
can be expressed as a summation, namely, 

E *=l E l B »,l exp(l^) I (43) 

K= 1 A 

and therefore the radar cross section is 



o = IE/To 7 exp(i^)|= 

K-l A 

r ,47 \d kw , . ,471 d K .. l7 

= l£ [cos( — *-*) + jsin( — 5 — ) ] I 2 

K - 1 



(44) 



where 

a*, is the crossing section of individual scattering 
elements , 

n the number of scattering elements, 
d K the range of the Kth element, and 
A, the wavelength. 

In order to relate the proceeding to a real target two 
assumptions shall be made: 
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(1) On a short term basis 4nd K /X is a random 
variable which can take on any value from 0 to 2k with equal 
probability , 

(2) The individual scattering elements have equal 
radar cross sections, that is, ak=c 0 . 

Therefore the random variables can be expressed as the x 
component (a 0 ) <1/2| cos (4Kd k /’K) and the y component (G 0 ) 1/2 
cos (4 Kd k /\) , the problem of determining the probability 
density function of <7 is identical to the problem of 
determining the distance moved from an origin in a two 
dimension walk problem of n steps, where the length of each 
step is (G 0 ) 1/2 and the direction of each step is perfectly 
random . 

The probability of having a component (x,y) after n 
steps where n is a comparativity large number is given by 



w(x,y) dx dy 



exp [-U 2 +y 2 )/na 0 ] ^ 

Kno 0 



( 45 ) 



Converting from rectangular coordinates to polar coordinates 
yields 

w(i?,0) dR dd = — exp[-(i? 2 ) /no Q ] dR dQ ( 46 ) 

Kno 0 

The marginal distribution of R, obtained by integrating the 
preceeding with respect to 0,is 

w(p) dR = —exp[-{R 2 )/na Q ] dR ( 47 ) 

no 0 
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Now, since C=R 2 =x 2 +y 2 and dO=2R dR Eq.(47) becomes 



w{o) = 



exp [-o/noJ 



no r 



for a > 0 



( 48 ) 



w(O)=0 for O<0. This demonstrates that the probability 
density function of the cross section is exponentially 
distributed with an average cross section <5=na 0 . Thus, the 
average total echoing area is the sum of the individual 
echoing areas of the individual elements. 



a. Case 1 Scan-to-scan fluctuations of exponential 
density function. 

The scan-to-scan fluctuations of an exponential nature can 
be applied to targets such as jet aircraft when a radar having 
a fairly high pulse repetition rate and scan rate is employed. 

For case 1, Swerling derived the probability of detection 
by initially defining a target model where the received signal 
power is exponentially distributed, namely, 

w( x ,x) M exp [ -x/x] toI x> 0 

^ ( 49 ) 

= 0 for x < 0 

where x = input signal-to-noise ratio. 

x = average of x over all targets fluctuations. 

The characteristic equation for the probability density 
function resulting from the integration of N pulse returns 
from an exponential fluctuating target, when there is complete 
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correlation between pulses, is derived from the non- 
fluctuating case (Swerling case 5) characteristic function as 
follows : 



< JBL) 

~-Nx'p+l‘ 

c N =[c x ) N = e 



(P+1 ) N 

C N (p) = / w(x,x) — — dx 

Jo P+1 w 



~ Nxl -pk ] 



(p+1)"' 1 [l+p{l+Nx)] 



( 50 ) 



For noise only, the characteristic equation is the same as 
for case 5 . 






(1+p) w 



( 51 ) 



The probability density function f n (y) obtained by using 
Campbell and Foster 431, p.44, is the same as the probability 
density function of the nonfluctuating case. i,e, 



f n (Y) = 



_y »- 1 e -y 



( 52 ) 



(1 +N) ! 

For signal plus noise, the probability density function 
f s + n (Y) was obtained by using Campbell and Foster pair 581.7, 
p.64. For N> 1, the density function f s . n (Y) is 
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/ • 1 
— (p+ l) w_1 ( 1 +p ( 1 +Nx) ] 



1 ^ -y/ ( i *//x) 

7® 



eJ'^py dp 



r (jv-i) ( - M e l ) 

i + wx 



- [r(w-i)-r(N-i / -i^] 
x \ W_1 1 + Mx 



N-2 



=-L e -y/n*m ( i +_±_ ) " ''(1-. 

Nx Nx 

= _L e -y/ (1 *^ (i+_l_) W 2 j( 



r(N-l,_^-y) 



Hn-d 



1 +NX ) 



y 



nx 



Nx 



l + _L yTFI 

Nx 



,N~ 2) 



(53) 



where T(v,z) is the incomplete gamma function and Ifu,pjis the 
incomplete Pearson gamma function. 

For N=l, the probability density function can be obtained 
from Campbell and Foster pair 438, p.45 as 

£.. a [Y)°eW*/{l+xl (54) 

therefore the probability of detection can be obtained by 
integrating Eq(54) from 0 to Y b and the result is given by 
Swerling as: 

1 ~P D =f* b f B . n (y) dv-I [ , N- 2 ] - ( 1 +n£) f ( Y b ) (55) 

From above, P D can be obtained as 



=1-I[— (N-2)) +(i + ^) W ‘ 1 






!(■ 



(1+ (l/Nx)’JTFI) 



Nx 

,N- 2 )e- r * /<1 * J ® 



(56) 
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Another approximation can be made by inspecting Eq(56) and 
let Nx>> 1 and P fa <<l such that the Pearson incomplete gamma 
functions are close to unity. Then Eq(56) can be rewritten 
approximately as 



1 

D n s (l + _ ) e 

Nx/ 2 



Nx>l 
Pfa< 1 



( 57 ) 



Taking the logarithm of Eq(57) and using a series expansion of 
each term, 



lnP d = (N-l ) In ( 1 +— ) £ 

NX {m d + J.) 

Nx 

e(N-l) [ 1 _ 1 1 - 1 1 

( Nx) 2 (Nx/ 2) 2 3 (Nx/ 2) 3 

i 

Nx/ 2 (Nx/ 2) 2 (Nx/ 2) 







( 58 ) 



Taking only the first terms and using Eq(40) yields 

In p D *-^(Y b ~N+l) 

** ( 59 ) 

+i] 

Nx 

b. Case 2 Pulse-to-pulse fluctuations of eicponential 
density function. 

Pulse-to-pulse fluctuations of an exponential nature apply 
to targets such as propeller-driven aircraft (if the 
propellers contribute a significant portions of the echo 
area) , to targets where small changes in orientation would 
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establish significant changes in echoing area (such as long 
thin subjected to a high-frequency signal), or to targets 
viewed by radars with sufficiently low repetition rates. 

For case two, the received signal return still belongs to 
exponential distribution. For this case the signal is 
completely decorrelated (pulse to pulse fluctuations exist) , 
the probability density function for noise only (signal- to- 
noise ratio equals to 0) still the same with Case 1. The 
probability of density function of signal plus noise for case 
2 can be obtained with the same procedures, but have the pulse 
to pulse integration. The characteristic equation for the 
single pulse is obtained by letting the characteristic 
equation of case 1 to be N= 1 



Ci(p) = 



( 60 ) 



[l+p(l+x) ] 

For N pulses which is completely decorrelated (independent), 
the characteristic equation is just the Nth power of single 
pulse 



C N ip) 



[Cjlp)]" = 



[ — 1 — ]» 
l+p(l +5c) 



( 61 ) 



Therefore the corresponding probability of density 
function obtained by using Campbell and Foster "Tables of 
Fourier transforms" pair 431, p.44 yield: 
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1 



( 62 ) 



= /-- l T^W ] e dy * u»a'w-i)' 



The probability of detection can be expressed as 



p D = 

= x- 1 [ Yb y_ N ' x c -v/n*xi dy_ 
[N-l) ! Jo i +x l +x 



= 1-J[- 



(1+3?) y/N 



, (N-l) ] 



( 63 ) 



Another approximation can be made by applying that the 
Gaussian family is closed under linear operation and by using 
Gram-Char lier series expansion (Appendix A) so that 



= (-j) 



dCy{p) 



d P 



lp=Q 



= N(l+x) 



d 2 Cy{ P ) 



m 2 = ( -J ) 2 lp-o = W+D (l + *) 2 - 

dp w ^ 

Variance : o 2 = -ml =N{ 1 +x) 2 



( 64 ) 



The approximate probability density function then is 



f B +n(Y)=- 



•/TkN{1+5c) 



-[r-Ar(i+x )] 2 
e 2N(l*"x) 2 



fjY) =_l_e<-™> 2 /2w 



( 65 ) 



Therefore the probability of detection and false alarm can 
be obtained by integrating from the threshold level y b to °° 
and are represented as: 
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( 66 ) 



Po=<t>[ 



(Y b ~N( 1+3?) ) 
v/3v(l+3c) 



^a=«M 



s/N 



2. SWERLING'S CASE 3 AND 4 

For Swerling's case 3 and case 4, the radar cross section 
can be described by chi-square distributions with four degrees 
of freedom. The density function is commonly associated with 
tabilized missile tankage and can be expressed as: 

A -(—) 

w{x,x) = — e * x>0 (67) 

x 2 



a. Case 3 Scan-to-scan fluctuations with a chi-square 
density function with four degrees of freedom. 

The characteristic equation for case 3 which represents 
the condition of complete correlation (scan-to-scan, no pulse- 
to-pulse fluctuation) is given by : 



C(p) ■ /• 



4x 

p 



-< 4 ?) 



(P+1)" 



dx 



(p+l) 



2-N 



[iTd + fn 



( 68 ) 



The characteristic equation and the probability of density 
function for noise only remain the same as in case 1 and 2 



C N (p) = 



(p+1) N 



e j2npy dp = 



y"- 1 e' y 



(69) 



(p+1)"' ' (N- 1) ! 

For signal plus noise the probability of density function 
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can be expressed as : 



f 8 .n(y) = /" (p * 1 ) 2 i" a e j2npy dp 

J -® r - . _ / - . A7n \ i 



[i +P (i + ito)] 2 

y *- 1 e~ y 

(AT-1) ! [(i+(JVx/2)] 



(70) 



(2 , N, 



1+ (2/Nx) 



where the density function comes from Fourier transform pair 
581.1 in the Campbell and Foster tables given as 



ya-1 £-(p o) Yf 2 x 



(s+p) “ +v (s+o) “~ v r(2a)(p-o)“ 

^.-i/ 2 [(p-o) r] ,y< o 



(71) 



where 



[i+l - — *i 

i^U) = Z 2 e 2 F l (p -v +-i- , 2|l +1 ; Z) 



(72) 



Eq(72) can be simplified using the relationships for 
confluent hypergeometric function, that is 

F 1 (2, N; Z) = (Z+2-N) F 1 {1,N;Z) +N-1 

FAl.N-.Z) = e z Z' w+1 (N-l) ! J[ ,iV-2] 

vT^T-TT 



the 



(73) 



Swerling uses two identities relating to confluent 
hypergeometric functions in order to expand the proceeding 
into more familiar Pearson's form of the incomplete gamma 
function, I(u,p) which is defined in Eq (30), thus Eq (73) 
becomes 
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(!/•)= [l + (2/^x)J \_ N 2 y rr y n-'* 1 e (-y/i + Wx/ 2 )) 

(y) [ 1 + ( nx/ 2) ] 2 I[ [i +2 TUmu^ 



- (^~ 2) fit (2/ (flfr) l"' 1 j[ Z , N-2 ] ( 74 ) 

[ 1 + ( Nx/ 2) ] 2 [1+ (2/Nx) ] v^Fl 

Q-y! (i+Nx) + y W-1 e~ y 

[1+ (JVx/2) ] (W-l) ! 



The probability of detection can be obtained by 
integrating the density function from the threshold Y b to °» 
as : 



P D = 



y%' 2 e~ Yb 2 Y, 



N-2 



_ + y e - y »Il + ( m±z ) ”-* e Nx*2 

{N-2) ! Nx+ 2 fa m! Nx 



tf -2 



~2Y h 



7Y m 

N-2 -2 Y b ( ^ ) 



(75) 



[ 1 - l (N ~ 2 .).+ 2 J b ] [x-y e~2^Nx 2 ±NK ) ] ; 

Nx Nx+2 



Ni.2 



m-Q 



b. Case 4 pulse-to-pulse fluctuations with a chi- 
square density function with four degrees of 

freedom 

With no correlation, the characteristic equation for a 
single hit is obtained by letting N= 1 in Eq(52) 



C(p) 



1+P 

tl+p(H--|)] 2 



( 76 ) 



The characteristic equation for the sum of N pulse can be 
obtained by the Nth power of a single hit 
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(1+p) " 
l+p(l+x/2) ] 



(77) 



C w (p)=[C 1 (p)) w = 






The inverse Fourier transform of above equation is 
obtained from Campbell and Foster transform Pair 581.1,p.64, 
this yields the following expression for the probability of 
density function in the condition of signal plus noise 



f (y) = f“i* 

Wy; Jr* 3? 



-(4?) 



(1+p) " 



[l+p(l+x/2) ] 



2AT 



-y 



(78) 



yN 1 Q l + x/2 



(1+3?) 



2AT 



(N- 1 ) ! F 1 (-N,Ni X i 2 y) 
l+x/2 



The confluent hyper geometric function F 1 (-N,N;a) can be 
expanded as : 



F 1 (-N, N; a) 



= 1 + ZM.L 

N 1! 



(-N) (-W+1) a 2 „ (-AQ (-jy+1) (-A/+2) a 3 ^ 
W(W+1) 2! N(N+1) (N*2) 3! 



_ Y' (-1) k [NI/{N-K) ! ] a_f 
&S [(tf+tf-1) ! / (AT-1) ! ] tf! 



(79) 



From Eq(78) and Eq(79), the density function f s< . n (Y) can be 
rewritten as 



= y Af ~ 1 e 

^8+n 



- y 
l+x/ 2 //! 



W — 



y' / x/2 

( 1 +3?/ 2 ) 2W fe l+x/2 



* y* (80) 

[ (iV+iC-1) ! (JNT-/0 ! iC! ] 



The probability of detection is obtained by integrating Eq(80) 
from the desired threshold Y b to From the definition of 
incomplete gamma function, the probability of detection can be 
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written as: 



p d = fj f e +n (y) dy 

J Y b 

Y ( 81 ) 

w _ H - ■ : i N+K-l] 

= 1 _ Nl { x, K (l+x/2) y/N+K 

(l+x/2) w jfe 2 } K\ ( N-K)\ 

With Eq(40), the probability of detection can be approximated 
as : 



N 



Pn a 1- 



— ™ T (1) K 

[l+a/2) N bi 2 



_ I[^-h 



sfN (J)' 1 (P fa ) +N 



(l+g/2) JN+K 



,N+K-1] 



( 82 ) 



K\ ( N-K ) ! 



C. SEARCH RADAR DETECTION RANGE CALCULATION 

The Marcum- Swer ling theory represented by extensive sets 
of curves from the computer programs can be used to determine 
the detection range of a practical radar by introducing 
detection loss and others parameters. There are many types of 
detection losses which have been identified so far, and when 
these are considered, reasonable predictions of radar 
performance can be obtained. 

The radar's detection range can be determined by applying 
the desired signal-to-noise ratio determined from the Marcum- 
Swerling theory and the calculated detection loss, as given by 
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^nax - f 



( 83 ) 



P t G t G r X 2 o 

(4k ) 2 kTBF n L ( -^) 



] 1/4 



(m) 



where 

R^x = The maximum detection range in meters 
P t = Peak transmitter power in Watts 
X = Wavelength in meters 

G t r = Transmitter and receiver antenna power gains 
a = Average radar cross section in square meters 
k = Boltzmann's constant= 1.38 x 10* 23 J/deg. 

T = Effective system input noise temperature in degrees 
Kelvin (° K) 

B = Receiver bandwidth in Hertz . 

L = Detection system power loss factor 
F n = The receiver noise figure 

S/N= The smallest output signal-to-noise ratio 
When n pulse are integrated previous equation can be 
written as 



■^nax t ' 



P t G t G z X 2 on 
(4 %) 2 kTBF n L{^) n 



1/4 



(m) 



( 84 ) 



where the parameters are the same as that of Eq (83) except 
that (S/N) n is the signal-to-noise ratio of one of the n 
pulses that are integrated to produce the required probability 
of detection for a specified probability of false alarm. 
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III. SOFTWARE DEVELOPMENT 



This chapter describes the development of MATLAB programs 
for the efficient and accurate computation of probability of 
detection based on Marcum and Swerling theory of radar 
detection. The MATLAB source code is given in Appendix B, and 
complete programs are available from Professor G.S. Gill, Code 
EC/Gl, Naval Postgraduate School, Monterey, CA 93943. 

A. PROGRAM STRUCTURE 

The overall program structure is shown in Figure (5) . The 
structure is that of a main menu program which calls various 
submenu programs (mscurve.m number. m number. m ) as required. 
The submenu programs, when called, will then display the 
purpose of the subprograms. When called, the subprograms will 
call the function programs to do the actual computation. It is 
possible to exit the process from either the main program or 
from the subprograms or the function programs. The advantage 
of this format is that if the user wants to change one of the 
subprograms or function programs and wants to add other 
programs, this system can accommodate it. For each subprogram 
there is an error detect prevention and data entry double 
check function to inform the user and restart the process. 
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Figure 5 : Program structure 



1. User's guide and instruction 

To use these programs, the following MATLAB files have 
to be copied to the user subdirectory. A brief explanation of 
each is also given. 

•ms.m is the main menu program and gives brief descriptions of 
the M-S model and displays the main menu. 

•mscurve.m ,rmenu.m, and number. m are the submenu programs. 
These give the purpose of the programs when called. 
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•radar, m, radral.m, meyer.m, are the subprogram used to 
integrate swerll.m, swerl2.m, swrl3 .m, swerl4.m, swerl5.m 
function programs to calculate and plot the curves p d vs S/N, 
p d vs p fa and p d vs N, respectively. 

•rangel.m, range3 .m, range4.m, are the subprograms responsible 
for integrating and calculating the detection range. The user 
inputs the specified detection loss from detect. m. The 
rangel.m, range3.m, range4.m, will load the data from detec t.m 
and calculate automatically and display the detection curves 
on the screen. 

•point. m radthre.m, are the subprograms responsible for 
calculating and displaying the numerical results from the 
function program radpoint.m and marcum.m, respectively, 
•swerll, swerl2, swerl3, swerl4, swerl5, prob.m, thresh. m, 
threshm.m, bound. m, radpoint.m, noise. m, signal. m, are all 
function programs responsible for calculating the data from 
the subprogram and then return the numerical value. 

A 386 or 486 personal computer is suggested for greater 
speed. To start the program type ms and press [enter] at the 
MATLAB prompt . 

2. Printing graphical outputs 

Graphics output from all programs will be stored as meta 
files automatically. The operator can print out the desired 
graphics from MATLAB by typing Igpp <filename> [enter] . Once 
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the user restarts the programs, the previous meta files will 
be deleted automatically. 



3. Programs options 

Once the ms.m command is given, the first screen seen by 
the user is as shown in Figure 6 . 



%THIS PROGRAM IS DESIGNED TO ALLOW THE STUDENT TO 
%VARY THE PARAMETERS OF THE VARIOUS SWERLING CASES 
%IN ORDER TO STUDY THE EFFECTS. 

% 

% CASE #s DESCRIPTION 

% 1. Returned pulses are of a constant amplitude 
% over one scan, but are uncorrelated from 
% scan to scan. 

% 

% 2. Returned pulses are uncorrelated from pulse 
% to pulse and correlated from scan to scan. 

% 

% 3. Returned pulses are of a constant amplitude 
% over one scan, but are uncorrelated from scan 
% to scan. 

% 

% 4. Returned pulses are uncorrelated from pulse to 
% pulse and correlated from scan to scan. 

% 

% 5. The static case with constant S/N and pulse 

\ amplitude 



Figure 6 : Main menu descriptions 



This screen gives a brief description of the five target 
models. After pressing the "enter" key, the user will see a 
second screen as shown at Figure 7 . 
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MAIN MENU 

1) THE M-S CURVES 

2) NUMERICAL DETECTION PROBABILITY CALCULATION 

3) RANGE DETECTION CURVES 

Select a menu number: 

Figure 7 : Main menu 

At this point the operator has different choices to make 
depending upon what he/she wants to do. Once the operator 
chooses one of the items, the main menu program will transfer 
to the selected submenu in Figure 8. 



SUB MENU — [THE M-S CURVES ANALYSIS] 

1) PROBABILITY OF DETECTION vs. S/N 

2) PROBABILITY OF DETECTION vs. Pfa 

3) PROBABILITY OF DETECTION vs. N 

4) COMPARSION OF M-S CURVES 

Select a menu number: 

Figure 8: Submenu Screen 

From the submenu the user will have another set of 
choices. He can choose the item he wants to study. After he 
chooses one of the items, the following selected screen will 
be seen in Figure 9, Figure 10, Figure 11 and Figure 12. The 
user can follow the instructions on the screen to key in the 
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arguments he wants to study. Then a data input screen seen at 
Figure 13 will display the data for double check. After the 
completion of above procedure, a selected graphic output will 
appear on the screen as in Figure 14. After the graphics 
display, the user can press 'enter 1 to get the next screen as 
shown in Figure 15. At this point, the user can either choose 
to go back to the main menu or exit to print out the graphic 
display. A selected second submenu and third submenu are shown 
in Figure 16 and Figure 17; users can follow the same 
procedure to choose the items they want to study. 



\****** ******************** ******************************* 
% THIS PROGRAM RETURNS THE PLOTS FOR THE NUMBER 

% OF PULSES AND SWERLING CASE SPECIFIED IN THE PARAMETERS. 
% THE PLOTS WILL BE STORED IN METAFILES UNDER THE NAME 
% "RADAR. MET" FOR AN EASY PRINT OUT. 

% 

% (A) the swerling case number has to be determined now 

%★*★*****★******★****★★★★★**★★*★★★★★*★*****★★*★*****★*★★★★ 
echo off 

Enter the case number you want to study 



Figure 9: Subprogram Descriptions Screen (a) 



*********************************************************** 
(B) . The number of radar pulses the program is to 

integrate needs to be an integer between 1 and 600 
*********************************************************** 

cho off 

Number of Pulses to be integrated is n = 10 



Figure 10: Subprogram Descriptions Screen (b) 
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**********************************************************************^* 

(C) . The probability of false alarm rate curves (pfa) to be plotted 
must now be determined. Each choice of a pfa will result in a different 
curve being plotted on the graph. You need to choose the following; 

1. The smallest pfa curve to be plotted, pfamin = ? 

2. The largest pfa curve to be plotted, pfamax = ? 

3. The step size between pfamin aand pfamax, pfastep = ? 

If you wish to plot only one curve then enter the same value for 
pfamin and pfamax. 

The suggested default step size to use is that of PFASTEP = 10, 
which is quite sufficient. It is suggested that pfamin and pfamax 
be powers of 10 as that is the normal choice. 

★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★I 



Figure 11: Subprogram Descriptions Screen (c) 



************************************************************* 

The signal to noise ratio (S/M) in dB for which you wish to 
plot needs to be determined. The choices you need 
to make are; 

1. The smallest S/N point to be plotted, sdbmin « ? 

2. The largest S/N point to be plotted, sdbmax » ? 
Remember that S/N must be entered in dB. 

3. The stepsize between sdbmin and sdbmax «■ ? 

****************************************** i ***** ** 



Figure 12: Subprogram Descriptions Screen (d) 
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-6 

% CHECK YOUR PARAMETERS ! 

% 

echo off 

The case number you is 1.00 

The number of pulses you choice are 1.00 

The max false alarm probability you choice is 1.00e-12 

The min false alarm probability you choice arel.00e-12 

The pfa stepsize is 10.00 

The max signal-to-noise ratio you choice is 10.00 
The min signal-to-noise ratio -10.00 

The S/N stepsize is 1.00 

% 

%IF THE PARAMETERS ARE CORRECT PRESS 1 
%IF THE PARAMETERS ARE NOT CORRECT PRESS 2 
% 

Figure 13: Input Argument Checking Screen 



MARCUM-SWERLING CURVE 




Figure 14: Selected Result 
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If you want to go to the same submenu ENTER CHOICE = 1 
or 

If you want to go to the main menu ENTER CHOICE = 2 

or 

To exit this whole program PRESS RETURN 



Figure 15: Selected menu 



SUB MENU [ NUMERICAL CALCULATION] 

1) CALCULATION OF DETECTION PROBABILITY 

2) CALCULATION OF THRESHOLD LEVEL 



Figure 16: Selected second menu 



SUB MENU — [THE DETECTION RANGE ANALYSIS] 

1) RANGE DETECTION CURVES WITH DIFFERENT PROBABILITY OF FALSE ALARM 

2) RANGE DETECTION CURVES COMPARSION 

3) RANGE DETECTION CURVES USE THE DETECABILITY FACTOR 



Figure 17: Selected third menu 



B. Algorithms for M-S models 



The function programs in Figure 5 are based on the Marcum 
and Swerling detection theory. 

1. Function program algorithms 



a normal MATLAB function. The input arguments and output data 
are listed in Table [1] at the end of this chapter. All these 
programs were implemented with the help feature of the MATLAB 
environment. Typing "help < name of the function >" will 
explain how to use them independently. 



The detection programs start with the calculation of 
threshold y b by applying the equation (30) in Chapter 2. 



The false-alarm probability can be represented as a 
function of N and y b . To compute probability of detection for 
all the five cases Y b is required for a given p fa . 

Since Eq(85) is a finite power series, in order to find Y b 
for given P fal a recursive computation method has to be used 
rewriting Eq(85) 



The function programs can be executed independently as 



a. Threshold computation 




( 85 ) 
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( 86 ) 




=P fa (N~l,Y b ) +HN-1) 




=P fa (N, Y b ) +L(fl) 
=P fa (N,Y b ) +L(N- 1)-^ 



( 87 ) 



From Eq (86) and (87) 




( 88 ) 



For a single hit (N= 1) the false alarm probability and the 
relation with the following term (N=2) are 



Therefore this relation allow each term of the expansion 
to be based on the value of the proceeding term, therefore an 
algorithm can be formed to compute the values of the detection 
threshold y b and the number of integrated pulses N. 

The first procedure employed in the algorithm is to define 
an empirical threshold level Y 0/ then compute P(N,y 0 ) . On the 
basis of this empirically determined value y 0 , an empirical 
suggested Y b is given by 




L( 1) =Y b e' Yb 



( 89 ) 
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Y b =N-ftl+2 .3JZ{JL+JN-1) , where L = -logP fa 



( 90 ) 



This value is used as the starting point to compute 
P(N,y 0 ) then compare it against the desired value of P {a and 
the difference between Y N and W This value of correction Ay 
can be calculated by using Newton-Raphson method by noting hat 



ln-^ 



(N, Y n ) 



^i = V- 



fa 



e "Yjj / ( N-l) ! 
Pfa (X, Y N ) 



( 91 ) 



The procedure is repeated until the correction magnitude 
Ay/ (Ay+y) is indicated that Y b is within a sufficient 
accuracy. A computer independent algorithm notation is given 
on Figure 18. In Figure 18, the arrow implies a specification. 
The normal execution of statements is carried out line by 
line, starting at the top, but a branch may be designated by 
an arrow which results from the execution of a statement. A 
conditional branch is denoted by a colon statement, and the 
branch is executed if the comparison condition specified on 
the arrow is satisfied. Otherwise, the next statement in the 
sequence is executed. Notice in figure 18 that the program is 
terminated when the value Ay/ (Ay+y) is less or equal to 6. 
The value of £ can be assigned to be 10' e or 10‘ 12 . This 
accuracy should be sufficient for application. 
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- "logic (P/J 

n — y/H+2 . 3v/I(/E+/n-l) 

- 7o 

- exp ( -y) 



M «- 0 

M - M+l 

N : M 

Ym - Km* Y/Af 

Kms - Yms+Ym 

p *• Yms 

- (p/Ym) *ln(P/P f J 

- y+ay 



aY 

y 

c 



ill 



Exit Yk 



Figure 18: Algorithmic Program for Detection 
Thresholds, y b 



The MATLAB function for executing this process is named 
prob.m and thresh. m respectively. In Figure 18 £ is the 

smallest acceptable tolerance value. Ay/ (Ay+y) is compared 
with e, and when it is less than or equal to e, the 
computation will stop and return the value of false alarm 
probability. In thresh. m the smallest tolerance value was set 
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to be 10' 6 , this value is sufficient for desired accuracy. 



This recursive method will be used in all Swerling cases 
to find the threshold level y b 

h. Swerling case 1 algorithm 

Following equation is used to compute probability of 
detection for case 1. 



It is obvious to see that for N=1 the detection probability is 
e -(Yb / i*nx> ' a ig 0r ithm in computer independent notation which 
uses equation (92) to compute Pd: is shown on Figure 19 . Notice 
in Figure 18 that the detection threshold ( Y b ) is independent 
of the target fluctuation characteristics so that the 
algorithm given in Figure 18 is used to determine Y b for all 
target types. A MATLAB source code to compute probability of 
detection of case 1 targets is given in Appendix B. The name 
of this file is swerll.m. 



p d , =1-J[— W-2) 1 +(l + Az) N1 I( 

1 / AT -1 



v/JFl Nx 



[l+(l /Nx)fiFl] 



N- 2) e 




( 92 ) 
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Enter Y b ,N,X 




N * X 
Y b 
1 
0 

EXP (-Y) 

YM 

0 

M + 1 
H 

YM “ Y/M 
YMS + YM 
YMS 
A 
1 

P2 

Y/(l +1/X N ) 

(1 +1/X N ) N ' 1 

EXP (-Yb/ (1 + X M ) ) 

A1 * Aj o (1 - P 2 ) 

Pi + P, > Exit p Dl 

EXP [ -Yb/ (Xn + 1)] > Exit P Dl 



Figure 19s Algorithmic Program for Swerling 1 Target , P D1 



c. Swerling case 2 algorithm 

Eq(63) can be modified as: 
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( 93 ) 



A computer 
equation is 



Pn=l-I( — Z-—.N-1) 

* vma+xi 



^ (N ’TTk> 



independent algorithm to implement 
shown in Figure 20. 



the above 
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d. Swerling case 3 algorithm 



From Eq(75) , the equation for probability of detection for 
case 3 is as follows: 



P D3 e 



-- is ( 1 +-^y) 

2 ^ 






for N= 1 



r N-2 -Y h 



~2Y h 



Yre ' b + Pf {N-l , y J + { ) N ~ 2 e Ux+2 



(N-2)! N x +2 



{1 _ 2[NJ. ) + [l-p fa {N-l, ^1) ] ; NZ 2 



Nx 

YzNx 



Nx Nx+2 



Nx+2 



(94) 



An algorithm for above equation is presented in computer 
independent notation in Figure 21. 
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Enter Y b ,N,X 



XN c= 

Y «- 

— N 

A 4- 

YM 4 - 

YMS 4- 

M 4- 

M 4- 

— N-l 4- 

YM 4- 

YMS <- 

*-^P2 <- 

— 1 : 

A 4- 

P1 4- 

Y 4- 

^A1 4— 

A2 4- 

A3 4— 

P3 4- 

B 

n 2 

L 4- 

C 4- 

L 4— 

PO 

C 4- 

U B 4- 

P4 4- 

P 4- 

— G 4- 

H 4- 

P 4- 



N*X 

Y b 

1 

0 

EXP ( -Y) 1 

YM 

0 

M ♦ 1 — 

M 

YM Y/M 

YMS ♦ YM 

YMS 

A 

1 

P2 

Y (1 ♦ 2 / X H ) 

(1 ♦ 2/X*) m - j 
EXP ( -2Y b / (X*+2 ) ) 

1-2 (N-2) / (X* -4 2Y b / (X, 4 2) 

A1 * A2 * A3 * (1 - P2) 

1 

N 

N-2 

L 

L - 1 
L 

(C*L) 

C 

2Y b H ‘ l c^/BMXh^) 

PI ♦ P2 ♦ P3 > Exit P 0 , 

EXP -2Yb/ ( XN+2 ) 

Z XN Yb/ ( XN+2 ) 2 
G (1 -*• H) >Exit P 0J 



Figure 21: Algorithmic program for Swerling 3 Target, P D3 
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e. Swerling case 4 algorithm 

Swerling case 4 computer program modified expression can be 
obtained by applying the power series in Eq (30) into Eq(81) 



as 



Pn =1- 



m y 

(l+x/ 2 ) N fa> 



N [1 -P(N+k, l+x/2 ) 1 



-(• 



* 



ic! {N-k) ! 



( y 

^ NX ATI / \ * 

) 1 + _ /n-tf 2v! ( — ) 



1 + 

2 



1 + 



JJl-AT 



/n! 



K-Q 



Kl (N-K) 



(95) 



-(• 



Nx 



*E e 



( L_ ) m 

1 + ^ * JV! (-&* 

2 v- 2 j 



ml 



K- 0 



JC! (N-K) 



y 

«■! 1 ^’\ 

. . , e “ 2 y 2 

1+ M fa, m ' ^ 



2 



/r-o 



JC! (N-JC) ! 



( )« 

w . 1 + ^ 

e ^ J. ] 



m-2N 



ml 



2N-1 

-E 

m m Q 



, <— k>* 

‘T ^ 1 ltj T 



/n! 



2//-1 

-E 

m*N 



. r , ( Tjfc 1 " 

-73’ i*-f 



/n! 



f^T M ( 2 v K , l x, 

^ JC! (W-JO ! l+iVx/2 l+Nx/2 



NX 



‘PfM. 



l+Nx/2 



( — )* 

-< 7 ^’ 

)*£e ,T 



J 77*0 



/n! 



in-// 



ti-E 



itf! 



(• 



NX 

2 \r 



)*(■ 



^ JC! (JV-JC) ! l+Nx/2 l+Nx/2 
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)*-*] 



with this expression the detection probability computation for 
this case therefore can be simplified according to equation 
(95) and avoid the computation of infinite power series. A 
computer independent algorithm is at Figure 22. 
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f. Swerling case 5 algorithm 

In case 5, the expression of Marcum steady target model in 
Eq (38) can also be modified by substituting the infinite 
power series for the modified Bessel function I N (X) given by 




(~) 2K 

2 

K\ (N+K) ! 



( 96 ) 



and interchanging the order of summation and integration, one 
obtain the probability detection for case 5 as: 



Po5=e'" x ) 



K=0 



(Nx) * 
K\ 



N-1*K 



[ £ e- Y ±~] 
*-*- m! 



m=0 



( 97 ) 



By interchanging the order of the summation results in a 
efficient computation representation as: 



w-i 



/ 7?=0 



m 






rr.Jl 

ml 



m-N 



K= 0 



Nx (NX) 
K\ 



( 98 ) 



The right hand side of equation (98) has two terms, the latter 
term is a infinite summation of power series, therefore a 
recursive evaluation is necessary for computing Pvs- 

The method to handle this infinite power series is to 
separate above equation into two terms. Let 

Pds= p l + Po> 

where L represent the integrated pulses where L represent the 
integrated pulses . Let L k N, then P L can be represented as 
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( 100 ) 



N - 1 rr ni L y m m-N . \ k 

e -y b Y b a -r> Y » 1 1 



P L =V e - y »^+y e' ¥b — (l-y 

L ^ m ! m ! 4r*L 



/J7— 0 



m=N 



777 ! 



/r=o 



AT! 



the error term P„ can be represented as 



'.-5>- r £d-E 



m-N 



m=L+l 



17i ! 



/c=o 



,-wx [NX) K 
K\ 



if let P L =P D5 then the truncated error term will be 



*-Z ^<1-E 



K-N 



m»L+ 1 



/n! 



Jc=0 



(Nx) k e 
k ! 



k„-Nx 



m=0 K ‘ x=0 



J5c! 



Let L be large enough such that 



L+l-tf 



m=0 x=0 ^ ' 



Therefore 

An upper bound can be found to limit the desired 
and avoid the unnecessary computation. The 
independent algorithm is in Figure 23 . 



( 101 ) 



( 102 ) 



(103) 

(104) 

accuracy 

computer 
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2. Solving for round off error 

Round off error can take place during the calculation of 
case 5. That is because during the calculation of probability 
of detection, if the signal- to-noise ratio x is zero, the 
output detection probability will become the probability of 
false alarm, then the criterion in Eq (103) will become 



L 

I 

171 -Q 



-yY k , ,, L ^ N (Nx) k e' Nx s 



(L)* ( 1 - £ 



Kl 

L+l-N 



k-Q 



kl 



so 



k-Q 



(105) 



Since we can not set the minimum acceptable error to be 
zero, the computation process will not terminate. Due to this 
reason, trade-off between the accuracy and the maximum value 
of output data should be made. An algorithm to compute the 
threshold y b to deal with the round off error is in Figure 24. 
The corresponding MATLAB file name is Marcum. m which can be 
used to compute the maximum acceptable threshold level y b when 
the value of signal- to-noise ratio is zero and the number of 
pulses is large. 
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Enter P f ,,Y 0 



Y 

M 4 - 
YLX 4- 
YLOG 4- 
TSUM 4 — 

M 4- 
M 4 — 
SUNK : 
YLX 4— 
TSUM 4- 
> 4- 





|Y-YJ : 
Y » 

Y 4- 

M 4- 

YLX 4- 
YLOG 4- 
Y-YLX : 
M 4- 
SUNK 4- 
YLX 4- 
TUM 4- 
Y b 4- 



P ^ 

0 

0 

LOG ( Y ) 

EXP(-Y) 

M+l 
N-l 
LOG (M) 

YLX+YLOG-SUNK 
TSUM+EXP ( YLX-Y ) 

Y+LOG (XP rA *TSUM) *TSUM/ EXP (YLX-Y) 

e 

Yl > Yb 

Yl 

0 
0 

LOG ( Y ) 

P 

M+l 
LOG (N) 

YLX+YLOG-SUNK 
EXP (YLX-Y) 

Yl > Yb 



Figure 24: Algorithmic Program for Detection Thresholds, Y b 
with round off error preventation 
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3. Solving for underflow and overflow 

As mentioned before, the maximum acceptable number for 
MATLAB is nearly e‘ 709 ' 782712893384040999999 , the underflow and overflow 
problem arise due to the fact that the power series in M-S 
equations has the form 




If the calculation for this power series is over the 
maximum or below the minimum number that MATLAB can represent, 
then the MATLAB will return a string named NaN or empty matrix 
( [ ] ) and the whole computation process will be meaningless . 
If the computation needed large integrated pulse number and 
large threshold or signal-to-noise ratio, a special effort 
has to be made to retain the accuracy. One method is let the 
series be of the form 

k k 

e~ y =exp [-y+kln ( y ) In ( N) ] =e'P (107) 
ki w-i 

and let the value e'^ be compared to the smallest acceptable 
value e - 709 - 782712893384040999999 in MATLAB, of course e - 709 - 782712893384040999999 
is much smaller than the tolerance number € assumed in the 
program. If P is greater than 709.7821... then increase the 
number K in above equation until it is less than 709.7821... 
and then start the summation process. This method will save 
MATLAB computation of power series with a large number and 
avoid the under flow problems. An algorithm for using this 
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method to compute the detection probability is in Figure (25) , 
and the MATLAB source code is threshm.m and is given in 
Appendix B . 




of dtection, P as with underflow preventation 
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4. Input arguments and output data 

The table below summarizes all input and output data for 
the programs . 



PROGRAM 


INPUT DATA 


OUTPUT DATA 


THRESH. M 


Pfa N 


Y b 


PROB . M 


Y b , N 


Pfa 


SWERL1.M 


P f .»N, x 


^Dl 


SWERL2 .M 


Pfa/N, x 


p D2 


SWERL3 . M 


Pf a ^ N # x 


PD4 


SWERL4 .M 


Pfa/ N, X 


Pd4 


SWERL5 .M 


Pfa/N, X 


P D5 


THRESHM.M.m 


N, Pfa 


Y b 


MARCUM. M 


Y b , N , x 


Pol ° r P fa 


POINT. M 


N, x 


Pd (numerical) 


RADPOINT.M 


N, x 


Pd (output data) 


RADAR .M 


Pfa/N, X 


CHART S/N vs Pd 


RADAR 1 . M 


Pfa / N , X 


CHART Pfa vs Pd 


MEYER. M 


Pfa/N, X 


CHART N vs Pd 


RANGE. 1 


Pfa/N, X , 


HART R vs Pd 
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RANGE. 3 


Pfa' N, X , 


CHART R vs Pd 


Table 1 : Input arguments and output dat 


a 



Where 

N=Number of pulses to be integrated 
X =average signal-to -noise ratio 

Pd= probability of detection 
p fa =probability of false alarm 
R= the desired detection range 
ns=Swerling case number 
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IV. RESULTS 



The M-S detection probability curves can be obtain by 
applying the MATLAB programs in Chapter 3 . The results can let 
the user evaluate the detection probability detection 
performance easily. The detection range curves can illustrate 
relationship between the detection probability and the 
detection range by input the properly detection loss. 

A. PROBABILITY OF DETECTION CURVES 

As mentioned previously, the M-S model arise from the 
different fluctuations of target cross section. MATLAB 
programs can be used to plot probability of detection versus 
per pulse signal- to-noise ratio for given probability of false 
alarm for five cases. It can also be asked to determine per 
pulse signal- to-noise ratio required for given P D and P fa . 
This required signal-to-noise ratio can be used to compute 
maximum detection range from the radar range equation. 

Fig 26 and Fig 27 compare the five target models for a 
false alarm of 10' fi , and the number of integrated pulse as N=10 
and N=100 respectively. When the probability is larger than 
0.33, all four cases in which the target cross section is not 
constant requires greater signal-to-noise ratio than the 
constant cross section case. This increase in signal- to-noise- 
ratio will cause a reduction in detection range. Therefore if 
the characteristic of the target cross section are not 
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properly take into account, the actual performance of the 
radar might not measure up to the performance which is 
predicted from the constant target cross section. A greater 
signal- to-noise ratio is required when the fluctuations are 
correlated pulse to pulse (case 1 and case 3 ) than when the 
fluctuations are uncorrelated pulse to pulse (case 2 and case 
4) . Fig 27 also indicates that when the number of integrated 
pulses is large, the case 2 and case 4 will approach to the 
constant target case (case 5). 




Figure 26: Probability of detection curves for five target 
models 
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models 

Since both the false alarm rate P fa and detection 
probability P D are specified by the system requirement, the 
radar designer computes required signal-to-noise ratio from M- 
S curves and uses this to compute the maximum detection range. 
The greater the number of pulses integrated, the greater is 
resulting overall signal-to-noise ratio. This results in 
greater probability of detection, but it will require longer 
scan time. Figure 28 and Figure 29 give plots of P d VS SNR 
for selected P fa 's for case 1 for N equal to 100 and 500 
respectively . 
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Figure 28: Probability of detection curves for Swerling case 
1 




Figure 29: Probability of detection curves for Swerling case 

1 67 



Figure 30 and Figure 31 are the Swerling case 2 plots of 
P d Vs SNR with 10 and 100 pulses integrated respectively. From 
Figure 30, signal- to-noise ratio of 4.8 db approximately is 
required to yield a probability of detection of 0.8 with 10 
pulses integrated and probability of false alarm of 10' 6 . 




Figure 30: Probability of detection curves for Swerling case 
2 
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Figure 31: Probability of detection curves for Swerling case 
2 



Figure 32 and Figure 33 show the Swerling case 3 with 10 
and 100 pulses integrated respectively. 
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3 




Figure 33: Probability of detection curves for Swerling case 
3 
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Figure 34 and Figure 35 show Swerling case 4 pulse to 
pulse fluctuation with 10 and 100 pulses integrated 
respectively . 




Figure 34: Probability of detection curves for Swerling case 

4 




Figure 35: Probability of detection curves for Swerling case 
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Figure 36 and Figure 37 are Swerling case 5 steady-state 



target with 10 and 100 pulses integrated respectively. 




Figure 36: Probability of detection curves for Swerling case 
4 
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B. THE DETECTIOON RANGE CURVES. 



The ASR-9 radar is used as a sample to compute maximum 
detection range. ASR-9 radar is designed for surveillance of 
the terminal area around airports. A significant feature of 
this S-band radar is its use of moving target detection (MTD) 
processor which provides coherent integration among other 
features. The parameters of the ASR-9 radar are given in table 
below. Also tabulated is the empirically determined search 
detection range on aim 2 target for the various Swerling 
models with a P D =0.9 and P fa =10' 6 . 



ASR-9 PARAMETERS AND RANGE CALCULATION 



P e -12oo kW 



ranges 



£■2900 MHz 
T-1.05 us 
G «1 

G t -33.5 db 
G r -33.5 dB 
NF-5 dB 
L-12dB 
0-1.3° 
0-75°/sec 
PRF-1200 pps 
Bdopp" 150 Hz 
P n -0 .9 
P f% -10- 



R t -39.41 nm 
R, -52.41 nm. 
Rj-46.68 nm 
R t -57.29 nm 
R,-62.63 nm 



Losses 

Transmitter 

Receiver 

Mismatch 

Integrator 

Collapsing 

Beam Shape 

MTI 



2 dB 

2 dB 
1 dB 
1 dB 
1 dB 

3 dB 
.2 dB 



Total 12 dB 



Table 2 : ASR-9 PARAMETERS AND RANGE CALCULATION 
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models 




DETECTION RANGE IN NAUTICAL MILES 
Figure 39: Range detection curves 
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C. Numerical data output from probability of detection 
programs . 

A table of values below is included for the purpose of 
checking the programs . 



N 




casf 


x=3 . 162278 


x* 1 0 


x*31 .62278 


x* 100 


x*3 16 . 2278 


1 


13 . 8155055 


1 


.0361810926687 


.2848035870497 


.65475592794524 


.8721557722127 


.957383959769 

4 


2 


.0361810926687 


.2848035870497 


.6547559279452 


.8721557722127 


.957383959769 

4 


3 


.0202712698585 


.2918843939447 


.77944622707425 


. 9625961846032 


.995941497060 

5 


4 


.0202659173728 


.2918820913596 


.77944704886892 


.9652566475676 


.995941410232 

7 


5 


. 0045852017166 


.2480487019045 


.99722505883065 


. 9999995192496 


. 999999519246 
9 


3 


19.1291681 


1 


.1971745188135 


.5769062596383 


.83647025963833 


.9446918645347 


.982126063254 

8 


2 


.1630815180577 


.7468903012970 


.97820931643215 


.9990169462157 


.999965068946 

1 


3 


.1784365463374 


.6869964903258 


.94510656525289 


.9933271492948 


.999288753543 

2 


4 


.1368031960489 


.8340491394154 


.9977036152469 


.9999940489495 


.999999991764 

8 




5 


. 0888130400938 


.9727252372625 


.99999920934189 


. 9999992093418 


.999999209341 

8 


N 




cat* 


x*3 . 162278 


x* 10 


x»31. 62278 


x-100 


x* 316. 2278 


10 


32.7013405 


1 


. .4855434894507 


.7911151202538 


.9282416527967 


.9765960615308 


.992532970604 

2 


2 


.73389703346785 


.9989667533193 


. 99999988582414 


9999999999973 


1 


3 


.56937472704103 


.9139452495156 


.9800436879910 


.9988084719731 


.999877746289 

5 


4 


.78178938715194 


.9999629142286 


.99999999999768 


1 


1 


5 


.8355167894356 


.9999991524271 


.99999915242839 


.9999991524283 


.999999152428 

3 



N 




cat* 


x«3. 162278 


x*10 


x»31 . 62278 


x»100 


x»316 . 2278 


30 


63.5481801 


1 


.69852634183006 


.8917070566484 


.96429073448000 


.9885538115323 


.99639654382453 


2 


.99944729942719 


.9999999999989 


1 


1 


1 


3 


.8147880711125 


.9759019998584 


.99728588213879 


.99971801636709 


.99997145717056 


4 


.999931272763883 


1 


1 


1 


1 


5 


.999998744645914 


.9999992846098 


. 99999928460987 


.99999928460987 


. 99999284460987 
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V. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 



Radar detection theory is outlined in this thesis and 
computer programs are developed for probability of detection 
calculations in MATLAB . These programs are based on radar 
detection theory as developed by Marcum and Swerling. These 
programs give more accurate results than the commonly 
available programs based on the detectability method which is 
empirical is native. A user's guide and instruction are 
included in the thesis. Simple examples of the programs usage 
are illustrated in the thesis. 

Students can use the programs written for this thesis to 
investigate radar system performance . These programs are cost 
effective, convenient to use and easy to reproduce since they 
are run on personal computers that are readily available to 
students. In particular, the program can aid students by 
removing a major computational burden to allow him to perform 
real world detection calculations . These programs can also be 
used for calculations in radar development work. 



B . RECOMMENDATIONS 

Due to an overflow underflow problem the programs are 
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limited to integration of 600 pulses in the probability of 
detection calculations for all 5 cases. Most of the time this 
does not present any limitation as most radars integrate fewer 
than 600 pulses. However it may be desirable to remove this 
limitation by employing a Chernoff bound or other methods. 
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Appendix A 

The gram-Charlier Series is a series expansion of a 
probability density function in terms of the Gaussian 
function, its derivatives and the moments of the original 
density function. 

If P(x) is a probability density function that is nearly 
Gaussian and x has zero mean and unit variance, then the Gram- 
Charlier series expansion of P(x) is given by 



P(x) =V a i <J> (i) ( x ) where 

f=o 

. (A-l) 

4><°> ( x ) =<t>(x) = e~* 2/2 

s/Tii (A-l. a) 

4,(1) (y)- di »W 

dx 1 



Coefficient can be evaluated in terms of the Hermite 

polynomials : 

<l> i (x)=— ^ — e { ~* 2/2) H, (x) (A-2) 

SITU 7 

Some typical polynomial are 



H 0 {x) =1, 

# x {x) =x, 

H 2 {x) =x 2 -1, 
H 3 (x) =x 3 -x, 



{x) =x 4 =6x 2 +3 
H 5 {x) =x\-10x 2 +15x 
H 6 {x) =x 6 -15x i +45x 2 -15 
H n {x) =x 7 -21x 5 +105x 3 -105x 



A recursive relation for these polynomials is 



(A-3 ) 
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tf J+1 (X) =xH i {x) - iH i . 1 (x) 



(A-4) 



A particularly useful fact is that the hermite polynomials and 
the derivatives of the Gaussian function are biorthogonal . 
That is, 

f ty 1 (x) Hj (x) dx 

= (-l) i i! by (A-5) 

_ (-1) 1 i\ , i-j 
0 i*j 

To evaluate coefficient a it both side of Eq(A-l) are 
multiplied by Hi(x) and result is integrated from to 
With Eq(A-3),the result simplifies to 

a i = ( ~ 1) ^ f m p[x)H i {x) dx (A-6) 

The first three coefficients a 0 a x and a 2 can be found by 
substituting Eq(A-3)into Eq(A-6) 

a 0 = f P(x) dx=l 

J -00 

a 1 = -f xp(x)dx= 0 (A-7) 

* 2 =jrfl (x 2 -l)P(x)dX=i(l-l)=0 

With these results, the Gram-Charlier series expansion of p(x) 
in Eq(A-l) can be simplified to 

when the mean of random variable x is not zero and/or its 
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(A-8) 



p(x) =4> (x) +V (J) 1 (x) 

2^3 



variance o is not unity, the Gram-Char lier expansion of p(x) 
is obtained by expanding a related density function g(t) in a 
Gram-Charlier series, where 



(A-9) 

o 

By a simple transformation of variables it follows that 



p(x) 



sr(t) | 

Cfar \t-(x~x)/0 

dt 




x-X) 

o 



(A-10) 



If the Gram-charlier series expansion for g(t) is 



git) ( t) (A-ll) 

2~0 

then from Eq(A-lO) and Eq(A-ll), the expression for p(x) 
becomes 



P(x) =^;Y J 

° 2^0 



X-X j 

a 



(A-12) 



If both sides of Eq(A-12) are multiplied by Hj[{x- x )/ a] 
and the result integrated from -<* to °°, coefficient c L can be 
evaluated with relation (A-6) as 
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(A-13) 



Evaluating Eq(A-13) for the first few coefficients Ci yields 
00=1, C 4 = jt( a 4 " 3 ) 

c 1 =c 2 =0 , c 5 = --^-j-y a 5 ~10a 3 (A-14) 

c 3 = --j-j-a, c 6 = -^g-j-y (a 6 ~15a 4 +3 0) 



where is the ith central moment of p(x), normalized by ci : 

a i = —f [x-x ) i p(x)dx (A-14) 

- 1 a jl — 

The can be expressed in terms of the conventional moments 
of the distribution p(x); thus if m„ denotes the nth moment 
of p(x) that is 



m n = f x n p(x)dx 

J -00 



then os, through are related to noncentral moments 



_ m 2 -3m 2 m 1 +2ml 



a 3 = 

_ m i -im 2 m 1 +6m 2 ml-3ml 
a 4 - — 



«5 = 



_ /n 5 -5m 4 /n 1 +10;n 3 /n 1 2 -10/n 2 /n 1 2 +4/n 1 5 



«6 = 



_ m 6 - 6 / 775 / 7 ^ +1 5 /n 4 7 n 2 -2 0 m 2 ml +15 ir^ml -5ml 



(A-15 ) 

^ by 



(A-16 ) 



In summary, to obtain a Gram-Charlier series expansion, 
the moments 1 % are found from Eq(A-16) , or they can be 
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conveniently computed from the characteristic equation of the 
distribution. Next the 0^ are evaluated using Eq(A-17) . The 0^ 
are then substituted into Eq(A-14) to obtain coefficient Ci of 
Gram-Charlier expansion. 
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APPENDIX B 



%FILE NAME: MS.M 

clear 

clg 

clc 

echo on 



%THIS PROGRAM IS DESIGNED TO ALLOW THE STUDENT TO 
%VARY THE PARAMETERS OF THE VARIOUS SWERLING CASES 
%IN ORDER TO STUDY THE EFFECTS. 

% CASE # : DESCRIPTION 

% 1. Returned pulses are of a constant amplitude over one 
scan, but are uncorrelated from scan to scan. 

% 2. Returned pulses are uncorrelated from pulse to pulse 
and correlated from scan to scan. 

% 3 . Returned pulses are of a constant amplitude over one 
scan, but are uncorrelated from scan to scan. 

% 4 . Returned pulses are uncorrelated from pulse to pulse 
and correlated from scan to scan. 

% 5. The static case with constant S/N and pulse amplitude 

echo off 

pause 

clc 

k=menu ( ' MAIN MENU ' , ' THE M-S CURVES 1 , . . 

'NUMERICAL DETECTION PROBABILITY CALCULATION' , . . 

'RANGE DETECTION CURVES') 

if k==l 

ms curve 

elseif k==2 

number 

elseif k==3 

rmenu 

elseif k==4 

density 

end 

% Now restart the process with some different choices 
% Clear the workspace of unnecessary data to avoid 
conflicts 

clear i; clear sav; clear sdb; clear w, clear pfa, clear 
p; clear ns; 

clear i; clear savl; clear sdbl; clear w, clear pfal, clear 



p; 

% Check and see if any info is to change 
clc 

echo on 

% If you want to restart the process ENTER CHOICE = 1 

% 

% or 
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% 

% To exit this programPRESS RETURN 
echo off 

choice = input ( 'CHOICE = ' ) 
if choice==l 
ms 

else 

end 

clc 

echo on 

% You have chosen to exit the program 
% 

% Bye ! 
echo off 



% FILE NAME: RADAR. M 
clc 

format long 
echo 

^fr****************************************-**-**-***-********** 

% THIS PROGRAM RETURNS THE PLOTS FOR THE NUMBER 

% OF PULSES AND SWERLING CASE SPECIFIED IN THE PARAMETERS. 

% THE PLOTS WILL BE STORED IN METAFILES UNDER THE NAME 
% "RADAR. MET" FOR AN EASY PRINT OUT. 

% 

% (A) the swerling case number has to be determined now 

^************'**************'*****'*****'*'******-**W***-***W*-**** 

echo off 

ns= input ('Enter the case number you want to study '); 
echo on 

^fr*****************************************'**********-*-****'*-* 

★ 

% (B) . The number of radar pulses the program is to 

% integrate needs to be an integer between 1 and 

600 . 

echo off 

n=input (' Number of Pulses to be integrated is n = '); 
clc 

echo on 

^****«***********************************«.*.****'i r ************ 

% (C) . The probability of false alarm rate curves (pfa) 

to be plotted 

% must now be determined. Each choice of a pfa will result 
in a different 

% curve being plotted on the graph. You need to choose the 
following; 
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% 1. The smallest pfa curve to be plotted, pfamin = ? 

% 

% 2. The largest pfa curve to be plotted, pfamax = ? 

% 

% 3. The step size between pfamin aand pfamax, pfastep = 

7 

% 

% If you wish to plot only one curve then enter the same 
value for 

% pfamin and pfamax. 

% 

% The suggested default step size to use is that of 

PFASTEP = 10, 

% which is quite sufficient. It is suggested that 

pfamin and pfamax 

% be powers of 10 as that is the normal choice, 

echo off 

pfamin = input ( 'pfamin = ’) 
pfamax = input (' pfamax = ')• 
pfastep=input ( 'pfastep = ') 
clc 

echo on 

% ‘ (D) . The signal to noise ratio (S/N) in dB for which 
you wish to' 

% plot needs to be determined. The choices you need 

% to. make are; 

% 

% 1. The smallest S/N point to be plotted, sdbmin = ? 

% 

% 2. The largest S/N point to be plotted, sdbmax = ? 

% 

% Remember that S/N must be entered in dB. 

% 

% 3. The stepsize between sdbmin and sdbmax = ? 

echo off 

sdbmin =input('In dB the sdbmin is ') 
sdbmax = input ( ' In db the sdbmax is ') 
sdbstep= input ('In db the sdbstep is ') 
clg 
clc 

echo on 

% 

% CHECK YOUR PARAMETERS ! 



echo off 

fprintf('The case number you is %8 . 2f \n ' , ns ) ; 
fprintf('The number of pulses you choice are %8.2f\n',n); 
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fprintf('The max false alarm probability you choice is 
%8 . 2e\n ' , pfamax) ; 

fprintf('The min false alarm probability you choice 
are%8 . 2e\n ' , pfamin) ; 

fprintff'The pfa stepsize is %8 .2f\n' ,pfastep' ) ; 
fprintf('The max signal- to-noise ratio you choice is 
% 8 . 2 f \ n ' , sdbmax ) ; 

fprintf('The min signal- to-noise ratio %8 . 2f\n sdbmin ) ; 
fprintf('The S/N stepsize is %8 . 2f\n sdbstep ') ; 
echo on 

% 

%IF THE PARAMETERS ARE CORRECT PRESS 1 
%IF THE PARAMETERS ARE NOT CORRECT PRESS 2 

% 

echo off 

choice=input ( 'CHOICE= ' ) ; 

if choice==2 

radar 

elseif choice==l 
end 

echo on 

% Now run the calculations for the above data and plot it 
echo off 

%delete radar. met 

format long 

nmax = 700; 

axis ( ' square '•) ; 

w= [sdbmin sdbmax 0 100]; 

axis ( w) ; 

if n < 1 | n > 600 | n-=fix(n), 

error (' Number of pulses input exceeds dimensions') 
end 

if ns < 1 | ns > 5 | ns ~= fix (ns), 

error (' Swerling case input does not correspond to 
allowable choices ' ) 
end 

pf a=pf amax; 

sdb = sdbmin : sdbstep : sdbmax ; 
kk=length (sdb) ; 
sav = 10 . ~ (sdb/ 10) ; 
if ns == 1 

% swerlingl 

while pfa >= pfamin 
for i= l:kk 

%i = 1 : sdbmax-sdbmin+1 
p ( i ) =swerll (pf a, n, sav ( i ) ) ; 
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p(i) =100*p(i) ; 

end 

plot (sdb, p, ' - ' ) ; grid; title ( ' SWERLING 1 ' ) ; 
xlabel ( ' (S/N) 1, signal- to-noise ratio, DB ' ) , 
ylabel (' Probability of detection in percent ') 
text(sdbmin + l/20*sdbmax, 95 ,['n = ' , int2str (n) ] ) 
text(sdbmin + l/20*sdbmax, 90,[*pfa = 

' , num2s tr (pf amax) , ' to ' , num2str (pfamin) ] ) 

text(sdbmin + l/20*sdbmax, 85,['pfastep = 

1 , num2str (pf astep) ] ) ; 

gtext( [num2str (pf a) ] ) ; 
hold on 

pfa = pfa/pfastep; 
end 
pause 

met a radar 
elseif ns==2 
%swerling2 
pfa=pfamax; 
while pfa > = pfamin 
for i= 1 : kk 

% for i = 1 : sdbmax-sdbmin+1 
p ( i) =swerl2 (pfa, n, sav( i ) ) ; 
p ( i ) =100*p(i) ; 
end 

plot ( sdb, p, title ( ’SWERLING 2’); 
grid, xlabel (' (S/N) 1, signal-to-noise ratio, DB'), 
ylabel (' Probability of detection in percent') 
text(sdbmin + l/20*sdbmax, 95 ,['n = ' , int2str (n) ] ) 
text(sdbmin + l/20*sdbmax, 90, [ 'pfa = ' , num2str (pfcuricix) , ' 
to ', num2str (pfamin) ] ) 

text(sdbmin + l/20*sdbmax, 85, [ 'pf astep = 

' , num2str (pf astep) ] ) 

gtext ( [num2str (pfa) ] ) ; 
hold on 

pfa = pfa/pfastep; 
end 
pause ; 
meta radar 
elseif ns==3 
%swerling3 
if n == 1 

error (' algorithm for swerling 3 does not work for n=l ' ) 
end 

while pfa >= pfamin 
for i=l:kk 

%f or i = 1 : sdbmax-sdbmin+1 
p( i) =swerl3 (pfa,n, sav(i) ) ; 
p(i) =100*p ( i) ; 
end 
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plot ( sdb, p, ' - ' ) , title ( ' SWERLING 3 ' ) ; 
grid, xlabel ( 1 ( S/N) 1 , signal- to-noise ratio, DB ' ) , 
ylabel (' Probability of detection in percent') 
text(sdbmin + l/20*sdbmax, 95 ,['n = ' , int2str (n) ] ) 
text(sdbmin + l/20*sdbmax, 90 , ( ’pfa = ' , num2s tr (pf amax) 
to ' , num2str (pf amin) ] ) 

text(sdbmin + l/20*sdbmax, 85 , [ ' pf astep = 

' , num2str (pf astep) ] ) 

gtext ( [num2str (pfa) ] ) ; 
hold on 

pfa = pfa/pfastep; 
end 
pause ; 
met a radar 
elseif ns==4 
%swerling4 

while pfa >= pfamin 
for i=l : kk 

%for i - 1 : sdbmax-sdbmin+1 
p(i) =swerl4 (pfa,n, sav(i) ) ; 
p ( i ) =100*p ( i ) ; 
end 

plot ( sdb, p, •-' ) ; title ( 'SWERLING 4 ' ) ; 
grid; xlabel ('( S/N) 1 , signal- to-noise ratio, DB ' ) ; 
ylabel (’ Probability of detection in percent') 
text(sdbmin + l/20*sdbmax, 95 ,['n = ' , int2str (n) ] ) 
text(sdbmin + l/20*sdbmax, 90, [ 'pfa = ' , num2str (pf amax) 
to ', num2str (pfamin) ] ) 

text(sdbmin + l/20*sdbmax, 85, [ 'pf astep = 

' , num2str (pf astep) ] ) 

gtext ( [num2str (pfa) ] ) ; 
hold on 

pfa = pfa/pfastep; 
end 
pause 

meta radar 
else 

%swerling5 

while pfa >= pfamin 
for i=l:kk 

% for i = 1 : sdbmax-sdbmin+1 
p(i) =swerl5 (pfa,n, sav(i) ) ; 
p(i) =100*p ( i ) ; 
end 

plot ( sdb, p, '-' ) , title ( 'SWERLING 5' ) ; 
grid, xlabel ('( S/N) 1 , signal-to-noise ratio, DB'), 
ylabel (' Probability of detection in percent') 
text ( sdbmin + l/20*sdbmax, 95 ,['n = ' , int2str (n) ] ) 
text(sdbmin + l/20*sdbmax, 90, [ 'pfa = ' , num2str (pf amax) 
to ', num2str (pfamin) ] ) 

text(sdbmin + l/20*sdbmax, 85 , [ 'pfastep = 
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■ , num2str (pf astep ) ] ) 

gtext ( [num2str (pfa) ] ) ; 
hold on 

pfa = pfa/pfastep; 
end 
pause 

meta radar 
end 

hold off 



%FILE NAME: SWERLl.M 
function pd=swerll (pfa, n, sbar ) 

% RETURNS THE VALUE OF THE PROBABILITY OF DETECTION FOR THE 
SWERLING1 

% CASE, GIVEN THE NUMBER OF PULSES n, THE PROBABILITY OF 
FALSE ALARM 

% PFA AND THE AVERAGE SIGNAL TO NOISE RATIO S/N. 

% EXAMPLE: SWERL1 ( PFA, N, SBAR) ; 

format long 
yb=thresh (pfa, n) ; 
y=yb/ (l+n*sbar) ; 
cte= ( 1+1/n/sbar ) ; 

pd=prob (n-1 , yb) +cte A (n-1) *exp ( -y ) * ( 1-prob (n-1 , yb/cte) ) ; 

return 

end 

% FILE NAME: SWERL2.M 
% 

***★★•*•**★★★★***★*★**************★*******★★*★****★★**•*•★**♦★★* 

% RETURNS THE VALUE OF THE PROBABILITY OF DETECTION FOR THE 
SWERLING2 

% CASE, GIVEN THE NUMBER OF PULSES n, THE PROBABILITY OF 
FALSE ALARM 

% PFA AND THE AVERAGE SIGNAL TO NOISE RATIO S/N. 

% EXAMPLE : SWERL2 ( PFA , N , SBAR ) ; 

^j**********************************************-*-****-*-******* 

function pd=swerl2 (pfa, n, sbar ) 

yb=thresh(pfa,n) ; 

y=yb/ ( 1+sbar ) ; 

pd=prob (n, y ) ; 

return 



% FILE NAME: SWERL3 .M 

% + + + + + + + + + + + '* + •* + + + + + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + + + '+ + + + + + + + + + 
% RETURNS THE VALUE OF THE PROBABILITY OF DETECTION FOR THE 
SWERLING3 

% CASE, GIVEN THE NUMBER OF PULSES n, THE PROBABILITY OF 
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FALSE ALARM 

% PFA AND THE AVERAGE SIGNAL TO NOISE RATIO S/N. 

% EXAMPLE: SWERL3 ( PFA, N, SBAR) ; 

function pd=swerl3 (pfa, n, sbar ) 
yb=thresh (pf a, n) ; 

y=i; 

for k=l:l:n-2; 
y=yb/k*y ; 
end; 

yl=2*yb/ (n*sbar+2) ; 
y=y*exp ( -yb) *yl; 

y2= ( (n*sbar+2 ) /n/sbar) ~ (n-2 ) *exp ( -yl) * ( 1-2* (n-2 ) /n/sbar+yl ) ; 
y2=y2* ( 1-prob (n-1 , yb*n*sbar/ (n*sbar+2 ) ) ) ; 
pd=y+pr ob ( n- 1 , yb ) +y2 ; 
return 

% FILE NAME: SWERL4.M 

^■fr********************************************************** 

% RETURNS THE VALUE OF THE PROBABILITY OF DETECTION FOR THE 
SWERLING4 

% CASE, GIVEN THE NUMBER OF PULSES n, THE PROBABILITY OF 
FALSE ALARM 

% PFA AND THE AVERAGE SIGNAL TO NOISE RATIO S/N. 

% EXAMPLE: SWERL4 ( PFA, N, SBAR) ; 

function pd=swerl4 (pfa, n, sbar ) 
yb=thresh (pfa, n) ; 
v=2 *yb/ ( sbar-t-2 ) ; 

ZK= ( 2 / ( sbar+2 ) ) ~n; 

YM=exp ( -v) ; 

YMS=YM; 
for M=1 : n-1 
YM=YM*v/M; 

YMS=YMS + YM; 
end 

SUM=YMS ; 

ZKS=ZK ; 

for M=n: (2*n-l) 

YM=YM*v/M; 

SUM = SUM + YM* ( 1-ZKS ) ; 

ZK = sbar*ZK* (2*n-M) / (2* (M-n-t-1) ) ; 

ZKS = ZK S + ZK; 
end 

pd=SUM; 

return; 

%FILE NAME : SWERL5 .M 

% RETURNS THE VALUE OF THE PROBABILITY OF DETECTION FOR THE 
SWERLING5 
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% CASE, GIVEN THE NUMBER OF PULSES n, THE PROBABILITY OF 
FALSE ALARM 

% PFA AND THE AVERAGE SIGNAL TO NOISE RATIO S/N. 

% EXAMPLE : SWERL5 ( PFA , N , SBAR ) ; 

function pd=swerl5 (pf a , n, sbar ) 
yb= thresh (pf a , n) ; 
z=n*sbar ; 

YM=exp ( -yb) ; 

OK = 0; 
k=n ; 

YMS=YM; 
for i=l : k-1 
YM=YM*yb/i; 

YMS = YMS + YM; 
end 

SUM = YMS ; 

XB=exp ( -z ) ; 

XBS=XB ; 
while OK == 0 
YM=YM*yb/k ; 

YMS = YMS + YM; 

SUM= SUM+ YMM1-XBS); 
el=l-YMS ; 

XB =XB*z/ (k-n+1) ; 

XBS = XBS + XB; 
e2=l-XBS; • 
er=el*e2 ; 
if er <= le-6 
OK=l ; 
end 

k = k + 1; 
end 

pd = SUM; 
return; 

% FILE NAME: THRESH. M 

%* RETURNS THE THRESHOLD yb FOR A GIVEN PROBABILITY OF FALSE 
ALARM (pfa) 

%* AND A GIVEN NUMBER OF PULSES 
%* EXAMPLE: THRESH ( Pfa , N) . 

^★★***************W*****-****-**********<#*-******************-** 

function y=thresh (pf a, n) 

format long 

if n < 1 | n-=fix(n), 

error ('Wrong input number of pulses'); 
break 
end 

if n==l 

y=-log(pfa) ; 
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return 

end 

1 = -loglO (pfa) ; 
n2=sqrt (n) ; 12=sqrt ( 1 ) ; 
y=n-n2+2 .3*12* (12+n2-l) ; 
comp=le-6 ; 
ratio=l; 

while comp <= ratio 
p=prob (n, y ) ; 
ym=l; 

for k=l:l:n-l; 

ym=ym/k*y ; 
end; 

ym=ym*exp ( -y ) ; 

dely= (p/ym) *log (p/pfa) ; % correction magitude 

y=y+dely ; 

ratio=abs ( dely /y ) ; 

end 

return 

% FILE NAME: PROB.M 

%* THIS PROGRAM RETURNS THE FALSE ALARM PROBABILITY P(n,yb) 
as in * 

%* SUMMATORY (yb A k/ fact ( k) ) for k=0 to n-1, AND IS A COMMON 
FUNCTION 

%* FOR THE PROBABILITY OF DETECTION COMPUTATIONS 
★ 

%* EXAMPLE: PROB ( N , YB ) 

★ 

function p=prob(x,y) 
if x< 0 | x-=fix(x) , 

error ( 'Number of pulses should be integer and greater than 
zero ' ) , break, end 
inf ==1.7 97693 13 4862069 *10^308 ; 

E=709 .782712893384040999999; 

p=0 ; 
if x==0 
return; 
elseif x==l 
p = exp(-y) ; 
return; 
else 

% [if x> 1 ] 

t=l; 

for k=l:l:x-l; 
t=t/k*y; 

P=t+p; 

end 

p= (p+1) *exp ( -y) ; 
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end 

return ; 

%File name: THRESHM.H 
function yb=threshm(pfa, n) 
tol=le-6 ; 

E=709 .782712893384040999999; 

1 = -loglO(pfa); 

n2 = sqrt(n) ;12 = sqrt(l) ; 

y=n-n2+2 .3*12* ( 12+n2-l) ; 

plog=logl0 (pfa) ; 
xpfa= 1/pfa; 
ok=0 ; 

while ok==0 

if y < E 

m=0 ; 
ylx=0; 

ylog=log (y ) ; 
tsum=exp ( -y ) ; 

for m = 0:l:n-l 
m=rrw- 1 ; 

svink=log(m) ; 

ylx= ylx+ylog -sunk; 

tsum=tsum+exp (ylx-y ) ; 

end 

yl=y+log (xpf a*tsum) * tsum/ exp (ylx-y ) 

if abs(y-yl) < tol 
yb=yl ; 
ok=l ; 
return 

else 

y=yi ; 
ok=0 ; 

end 



else 

m=0 ; 
ylx=0 ; 

ylog=log(y) ; 
okk=0 ; 

while okk==0 



m=m+ 1 ; 

sunk= log (m) ; 
ylx=ylx+ylog-sunk ; 
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if y-ylx <= E 



tsum=exp (ylx-y ) ; 
okk=l ; 
else 

okk=0 ; 

end 

end 

if m >= n-1 

yl=y+log (xpfa*tsum) * tsum/exp ( -y+ylx) ; 
if abs(yl-y) < tol 
yb=yl ; 
ok=l ; 
break; 
return 

else 

y=yl ; 
ok=0 ; 

end 



else 

for m=m: 1 ; n-1 

sunk=log (m) ; 
ylx=ylx+ylog-sunk ; 
tsum=tsum+exp ( -y+ylx) ; 

end 

yl=y+log (xpf a*tsxam) *tsum/exp ( -y+ylx) ; 
if abs(yl-y) <= tol 
yb=yl ; 
ok=l ; 
break ; 
return 

else 

y=yl; 
ok=0 ; 

end 

end 

end 

end 

%FILE NAME: MARCUM. M 
function [p, ymo] =bound (n, y, x) 
format long 
err=le-6 ; 

E=709. 7827128933 84040999999 ; 

%★***♦* start the program *************************** 

lamdal= (n/ (2*y) ) ^2+ (n*x) /y; 

lamda=l- (n/ ( 2 *y ) ) -sqrt ( lamdal ) ; 

auxl= ( -lamda*y ) + (n*x*lamda) / ( 1-lamda) ; 
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aux2=exp (auxl) / (1-lamda) ~n; 
if aux2 <= err 
aux3=n* (x+1 ) ; 
if y >= aux3 
p=0 ; return 
else 

p=l ; return 

end 

else 

ib=n-l; 
k=0 ; 
x=n*x; 
if x <= E 

xk=exp ( -x) ; 
xks=xk; 
m=0 ; 

if y<= E 

ym=exp ( -y) ; 
yms=ym ; 

% . 



for m=l :n-l 

ym=ym*y/m; 

yms=yms+ym; 

if err>(l-yms) 
p=yms ; 
ymo=ym 
x=x/n; 

return 

end 

end 
m=n- 1 ; 
sum=yms ; 
ymo=ym 
ok=0 ; 

while ok==0 
k=k+l; 
m=m+l; 
ym=ym*y/m; 
yms=yms+ytn; 
sum=sum+ym* ( 1-xks ) ; 
xk=xk*x/k; 
xks=xks+xk; 
ans= ( 1-yms ) * ( 1-xks ) ; 
if ans<= err 
ok= 1 ; 

end 

end 
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p=sum; ymo=ym; 



else 

ymly=0; 
ylog=log(y) ; 
okk=0 ; 

while okk==0 
m=m+l ; 
xmf =log (m) ; 
ymly=ymly+ylog-xmf ; 

if (y-ymly) <= E 

ym=exp ( -y+ymly ) ; 
okk=l ; 
else 

okk=0 ; 

end 

end 

if m <= ib 
yms=ym; 



% 



for ml=m: n-1 

ym=ym*y/m; 
yms=yms+ym; 
if err>(l-yms) 
p=yms ; 
ymo=ym 
x=x/n; 
ook=l; 
return 
end 

end 

sum=yms ; 
ymo=ym 
ok=0 ; 

while ok==0; 
k=k+l ; 
m=m+l; 
ym=yTO*y/m; 
yms=yms+ym; 
sum=sum+ym* ( 1-xrs ) ; 
xk=xk*x/k; 
xrs=xrs+xk; 
ans= ( 1-yms ) * ( 1-xrs ) ; 
if ans > err 
ok=0 ; 

else 

p=sum; 

x=x/n; 

break; 

return; 
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ok=l ; 

end 

end 

p=sum; ymo=ym; 



else % for m> ib; 

ymo=0 ; return 
yms=ym; 
sum=0 ; 
kd=m-n; 
if kd > x 

p=0; 

x=x/n; 
break ; 
return; 
else 

okkk=0 ; 

while okkk==0; 
k=k+l ; 
kk=xk*x/k; 
xks=xks+xk; 

while (k-kd) < 0 
okkk=0 ; 

end 

end 

ooook=0 ; 
while ooook==0; 

sum=sum+ym* ( 1-xrs ) ; 
xk=xk*x/k; 
xrs=xrs+xk; 
ans= ( 1-yms ) * ( 1-xrs ) ; 
while ans < err 
p=sum; 
x=x/n; 
ooook=l ; 
break ; 
return 

end 
k=k+l ; 
m=m+ 1 ; 
ym=ym*y/m; 
yms=yms+ym; 
ooook=0 ; 

end 

end 

end 

end 

else 

xklx=0 ; 
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xlog=log (x) ; 
ooook=0 ; 
while oook==0 
k=k+l; 

xklx==xklx+xlog-xkf ; 
while (x-xklx) > E 
oook=0 ; 

end 

end 

xk=exp ( -x+xklx) ; 
ib=ib+k; 

% 

xks=xk; 
m=0 ; 

if y< = E 

ym=exp ( -y) ; 
yms=ym; 

% 

for m=l : n-1 

ym=ym*y/m; 
yms=yms+ym; • 

if err>(l-yms) 
p=yms; 
ymo=ym; 
x=x/n; 
return 
end 

end 

sum=yms ; 
ymo=ym; 
ok=0 ; 

while ok==0 ; 
k=k+l ; 
m=m+l ; 
ym=ym*y/m; 
yms=yms+ym; 
sum=sum+ym* ( 1-xrs ) ; 
xk=xk*x/k; 
xrs=xrs+xk; 
ans= ( 1-yms ) * ( 1-xrs ) ; 
if ans<= err 
ok=l ; 

end 

end 

p=sum; 



else 

ymly=0 ; 
ylog=log (y) ; 
okk=0 ; 

while okk==0 
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m=m+ 1 ; 
xmf=log(m) ; 
ymly=ymly+y log-xmf ; 

if (y-ymly) < = E 

ym=exp ( -y+ymly ) ; 
okk=l ; 
else 

okk=0 ; 

end 

end 

if m <= ib 
yms=ym; 



ook=0 ; 

while ook==0 
while m-=ib 
m=m+ 1 ; 
ym=ym*y/m; 
yms=yms+ym; 
if err> ( 1-yms ) 
p=yms ; 
ymo=ym; 
x=x/n; 
ook=l; 
break; 
return 
else 
ook=0 ; 
end 

end 

end 

sum=yms ; 
ymo =ym ; 
ok=0 ; 

while ok==0; 
k=k+l ; 
m=m+l ; 
ym=ym*y/m; 
yms =yms +ym ; 
sum=suxn+ym* ( 1-xrs ) ; 
xk=xk*x/k ; 
xrs=xrs+xk; 
ans= ( 1-yms ) * ( 1-xrs ) ; 
if ans > err 
ok=0 ; 

else 

p=sum; 
x=x/n; 
break ; 
return; 
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ok=l ; 



end 

end 



else % for m> ib; 

ymo=0 ; 
yms=ym; 
sum=0 ; 
kd=m-n ; 
if kd>x 

p=0 ; 

x=x/n; 

break; 

return; 

else 

okkk=0 ; 

while okkk==0; 
k=k+l ; 
kk=xk*x/k; 
xks=xks+xk ; 

while (k-kd) < 0 
okkk=0 ; 

end 

end 

ooook=0 ; 

while ooook==0; 

sum=sum+ym* ( 1-xrs ) ; 
xk=xk*x/k ; 
xrs=xrs+xk; 
ans= ( 1-yms ) * ( 1-xrs ) 
while ans < err 
p=sum; 
x=x/n; 
ooook=l ; 
break ; 
return 

end 
k=k+l ; 
m=m+l ; 
ym=ym*y/m; 
yms=yTns+ym; 
ooook=0 ; 

end 



end 



end 

end 



end 

end 
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